setwd("E:\\5hmc_file\\2_5hmc_yjp_bam\\ASM\\bayes_pvalue_beta0")
file=read.table("pvalue_sig_statis.csv",head=T,sep=",")
sel=which(file$num1==2&file$twins==1)
file=file[!(file$num1==2&file$twins==1),]
file=file[file$num1>1,]
id=file$unitID

group1=c("X2B_X1T","M8_M7","M6_M5","M2_M1","M48_M47","M50_M49","M28_M27","M30_M29","M26_M25","M35_M36","M18_M17","M20_M19","M22_M21","M40_M39")

i=7
fn1=paste0(group1[i],".bayes_p.txt")
file=read.table(fn1,head=T,sep = "\t")
file$unitID=paste(file$chrom,file$position,file$ref,file$var,sep=":")
file$id=paste(file$chrom,file$position,sep=":")
file=file[file$id %in% id,]
file=file[file$normal_bayes_pvalue<0.05|file$tumor_bayes_pvalue<0.05,]
file=data.frame(unitID=file$unitID,normal_bayes_beta0=file$normal_bayes_beta0,normal_bayes_pvalue=file$normal_bayes_pvalue,tumor_bayes_beta0=file$tumor_bayes_beta0,tumor_bayes_pvalue=file$tumor_bayes_pvalue)
file$normal_group="nosig"
file[file$normal_bayes_beta0>0&file$normal_bayes_pvalue<0.05,]$normal_group="up"
file[file$normal_bayes_beta0<0&file$normal_bayes_pvalue<0.05,]$normal_group="down"
file$tumor_group="nosig"
file[file$tumor_bayes_beta0>0&file$tumor_bayes_pvalue<0.05,]$tumor_group="up"
file[file$tumor_bayes_beta0<0&file$tumor_bayes_pvalue<0.05,]$tumor_group="down"
file$pattern=paste0("normal_",file$normal_group,"-tumor_",file$tumor_group)
names(file)=c("unitID",paste(group1[i],names(file)[2:5],sep="_"),names(file)[6:8])
normal_down_tumor_down=file[file$pattern=="normal_down-tumor_down",][,1:5]
one_is_down=file[file$pattern=="normal_down-tumor_nosig"|file$pattern=="normal_nosig-tumor_down",][,1:5]
onedown_oneup=file[file$pattern=="normal_down-tumor_up"|file$pattern=="normal_up-tumor_down",][,1:5]
one_is_up=file[file$pattern=="normal_nosig-tumor_up"|file$pattern=="normal_up-tumor_nosig",][,1:5]
normal_up_tumor_up=file[file$pattern=="normal_up-tumor_up",][,1:5]


for(i in 8:10){
  fn1=paste0(group1[i],".bayes_p.txt")
  file=read.table(fn1,head=T,sep = "\t")
  file$unitID=paste(file$chrom,file$position,file$ref,file$var,sep=":")
  file$id=paste(file$chrom,file$position,sep=":")
  file=file[file$id %in% id,]
  file=file[file$normal_bayes_pvalue<0.05|file$tumor_bayes_pvalue<0.05,]
  file=data.frame(unitID=file$unitID,normal_bayes_beta0=file$normal_bayes_beta0,normal_bayes_pvalue=file$normal_bayes_pvalue,tumor_bayes_beta0=file$tumor_bayes_beta0,tumor_bayes_pvalue=file$tumor_bayes_pvalue)
  file$normal_group="nosig"
  file[file$normal_bayes_beta0>0&file$normal_bayes_pvalue<0.05,]$normal_group="up"
  file[file$normal_bayes_beta0<0&file$normal_bayes_pvalue<0.05,]$normal_group="down"
  file$tumor_group="nosig"
  file[file$tumor_bayes_beta0>0&file$tumor_bayes_pvalue<0.05,]$tumor_group="up"
  file[file$tumor_bayes_beta0<0&file$tumor_bayes_pvalue<0.05,]$tumor_group="down"
  file$pattern=paste0("normal_",file$normal_group,"-tumor_",file$tumor_group)
  names(file)=c("unitID",paste(group1[i],names(file)[2:5],sep="_"),names(file)[6:8])
  normal_down_tumor_down=merge(normal_down_tumor_down,file[file$pattern=="normal_down-tumor_down",][,1:5],by="unitID")
  one_is_down=merge(one_is_down,file[file$pattern=="normal_down-tumor_nosig"|file$pattern=="normal_nosig-tumor_down",][,1:5],by="unitID")
  onedown_oneup=merge(onedown_oneup,file[file$pattern=="normal_down-tumor_up"|file$pattern=="normal_up-tumor_down",][,1:5],by="unitID")
  one_is_up=merge(one_is_up,file[file$pattern=="normal_nosig-tumor_up"|file$pattern=="normal_up-tumor_nosig",][,1:5],by="unitID")
  normal_up_tumor_up=merge(normal_up_tumor_up,file[file$pattern=="normal_up-tumor_up",][,1:5],by="unitID")
}

library(dplyr)
library(ggplot2)
normal_up_tumor_up$group="normal_up_tumor_up"
normal_down_tumor_down$group="normal_down_tumor_down"
one_is_down$group="one_is_down"
file=rbind(normal_up_tumor_up,normal_down_tumor_down,one_is_down)
j=ncol(file)-1
file$mean.normal.beta0=rowMeans(file[,seq(2,j,4)])
file$mean.tumor.beta0=rowMeans(file[,seq(4,j,4)])
std <- function(x) sd(x)/sqrt(length(x))###算标准误差
file$std.normal.beta0=apply(file[,seq(2,j,4)], 1, std)
file$std.tumor.beta0=apply(file[,seq(4,j,4)],1,std)
file$normal.beta0.upper=file$mean.normal.beta0+file$std.normal.beta0
file$nomal.beta0.lower=file$mean.normal.beta0-file$std.normal.beta0
file$tumor.beta0.upper=file$mean.tumor.beta0+file$std.tumor.beta0
file$tumor.beta0.lower=file$mean.tumor.beta0-file$std.tumor.beta0
anno=read.csv("../20201107/CC.intersect.anno.hg19_multianno.csv",head=T)
anno$unitID=paste(anno$Chr,anno$Start,anno$Ref,anno$Alt,sep=":")
file=merge(anno,file,by="unitID")

normaldata=select(file,unitID,mean.normal.beta0,normal.beta0.upper,nomal.beta0.lower,avsnp150,group)
normaldata=arrange(normaldata,group,mean.normal.beta0)
normaldata$num=1:dim(normaldata)[1]
orders=select(normaldata,unitID,num)

p1=ggplot(normaldata,aes(y=mean.normal.beta0,x=num))+scale_x_continuous(breaks = normaldata$num,labels = normaldata$avsnp150,position="top")+
  geom_point(aes(color=group),size=1.5)+geom_errorbar(aes(ymin=nomal.beta0.lower,ymax=normal.beta0.upper,color=group),width=0.5)+labs(x="",y="con β0")+
  geom_hline(yintercept = 0,color="red",linetype="dashed")+scale_color_manual(values = c('#708090','red','blue'))+coord_flip()+theme_classic(base_size = 15)+
  theme(panel.grid.minor = element_blank(),panel.border = element_blank(),axis.line =element_line(size = 0.5))+theme(legend.position="none")
p1
tumordata=select(file,unitID,mean.tumor.beta0,tumor.beta0.upper,tumor.beta0.lower,avsnp150,group)
tumordata$na=""
tumordata=merge(tumordata,orders,by="unitID")

p2=ggplot(tumordata,aes(y=mean.tumor.beta0,x=num))+scale_x_continuous(breaks = tumordata$num,labels = tumordata$na)+
  geom_point(aes(color=group),size=1.5)+geom_errorbar(aes(ymin=tumor.beta0.lower,ymax=tumor.beta0.upper,color=group),width=0.5)+
  labs(x="",y="case β0")+geom_hline(yintercept = 0,color="red",linetype="dashed")+scale_color_manual(values = c('#708090','red','blue'))+theme_classic(base_size = 15)+
  theme(panel.grid.minor = element_blank(),panel.border = element_blank(),axis.line = element_line(size = 0.5))+coord_flip()+theme(legend.position="none")
p2
layout <- matrix(c(1, 1, 1, 1, 2, 2, 2), nrow = 1)
multiplot(plotlist=list(p1,p2),layout=layout)

CCid=file$unitID
original.normaldata=normaldata
orders=data.frame(unitID=original.normaldata$unitID,original.normaldata[,5:7])

###看单发样本中的
i=1
fn1=paste0(group1[i],".bayes_p.txt")
file=read.table(fn1,head=T,sep = "\t")
file$unitID=paste(file$chrom,file$position,file$ref,file$var,sep=":")
file=data.frame(unitID=file$unitID,normal_bayes_beta0=file$normal_bayes_beta0,normal_bayes_pvalue=file$normal_bayes_pvalue,tumor_bayes_beta0=file$tumor_bayes_beta0,tumor_bayes_pvalue=file$tumor_bayes_pvalue)
rt=file[file$unitID %in% CCid,]
names(rt)=c("unitID",paste(group1[i],names(file)[2:5],sep="_"))

for (i in 2:6) {
  fn1=paste0(group1[i],".bayes_p.txt")
  file=read.table(fn1,head=T,sep = "\t")
  file$unitID=paste(file$chrom,file$position,file$ref,file$var,sep=":")
  file=data.frame(unitID=file$unitID,normal_bayes_beta0=file$normal_bayes_beta0,normal_bayes_pvalue=file$normal_bayes_pvalue,tumor_bayes_beta0=file$tumor_bayes_beta0,tumor_bayes_pvalue=file$tumor_bayes_pvalue)
  ftmp=file[file$unitID %in% CCid,]
  names(ftmp)=c("unitID",paste(group1[i],names(file)[2:5],sep="_"))
  rt=merge(rt,ftmp,by="unitID",all=T)
}
file=rt
j=ncol(file)-1
file$mean.normal.beta0=rowMeans(file[,seq(2,j,4)],na.rm = T)
file$mean.tumor.beta0=rowMeans(file[,seq(4,j,4)],na.rm = T)
std <- function(x) {
if(length(x[!is.na(x)])==1){return (0)
}else(return(sd(x,na.rm = T)/sqrt(length(x[!is.na(x)]))))
}###算标准误差
file$std.normal.beta0=apply(file[,seq(2,j,4)], 1, std)
file$std.tumor.beta0=apply(file[,seq(4,j,4)],1,std)
file$normal.beta0.upper=file$mean.normal.beta0+file$std.normal.beta0
file$nomal.beta0.lower=file$mean.normal.beta0-file$std.normal.beta0
file$tumor.beta0.upper=file$mean.tumor.beta0+file$std.tumor.beta0
file$tumor.beta0.lower=file$mean.tumor.beta0-file$std.tumor.beta0
#file=merge(anno,file,by="unitID")
file=merge(file,orders,by="unitID")


normaldata=select(file,unitID,mean.normal.beta0,normal.beta0.upper,nomal.beta0.lower,avsnp150,group,num)

p3=ggplot(normaldata,aes(y=mean.normal.beta0,x=num))+scale_x_continuous(breaks = normaldata$num,labels = normaldata$avsnp150,position="top")+
  geom_point(aes(color=group),size=1.5)+geom_errorbar(aes(ymin=nomal.beta0.lower,ymax=normal.beta0.upper,color=group),width=0.5)+labs(x="",y="con β0")+
  geom_hline(yintercept = 0,color="red",linetype="dashed")+scale_color_manual(values = c('#708090','red','blue'))+coord_flip()+theme_classic(base_size = 15)+
  theme(panel.grid.minor = element_blank(),panel.border = element_blank(),axis.line =element_line(size = 0.5))+theme(legend.position="none")
p3
tumordata=select(file,unitID,mean.tumor.beta0,tumor.beta0.upper,tumor.beta0.lower,avsnp150,group,num)
tumordata$na=""

p4=ggplot(tumordata,aes(y=mean.tumor.beta0,x=num))+scale_x_continuous(breaks = tumordata$num,labels = tumordata$na)+
    geom_point(aes(color=group),size=1.5)+geom_errorbar(aes(ymin=tumor.beta0.lower,ymax=tumor.beta0.upper,color=group),width=0.5)+
    labs(x="",y="case β0")+geom_hline(yintercept = 0,color="red",linetype="dashed")+scale_color_manual(values = c('#708090','red','blue'))+theme_classic(base_size = 15)+
    theme(panel.grid.minor = element_blank(),panel.border = element_blank(),axis.line = element_line(size = 0.5))+coord_flip()+theme(legend.position="none")
p4
layout <- matrix(c(1, 1, 1, 1, 2, 2, 2), nrow = 1)
multiplot(plotlist=list(p3,p4),layout=layout)

###看健康样本中的
i=11
fn1=paste0(group1[i],".bayes_p.txt")
file=read.table(fn1,head=T,sep = "\t")
file$unitID=paste(file$chrom,file$position,file$ref,file$var,sep=":")
file=data.frame(unitID=file$unitID,normal_bayes_beta0=file$normal_bayes_beta0,normal_bayes_pvalue=file$normal_bayes_pvalue,tumor_bayes_beta0=file$tumor_bayes_beta0,tumor_bayes_pvalue=file$tumor_bayes_pvalue)
rt=file[file$unitID %in% CCid,]
names(rt)=c("unitID",paste(group1[i],names(file)[2:5],sep="_"))

for (i in 12:14) {
  fn1=paste0(group1[i],".bayes_p.txt")
  file=read.table(fn1,head=T,sep = "\t")
  file$unitID=paste(file$chrom,file$position,file$ref,file$var,sep=":")
  file=data.frame(unitID=file$unitID,normal_bayes_beta0=file$normal_bayes_beta0,normal_bayes_pvalue=file$normal_bayes_pvalue,tumor_bayes_beta0=file$tumor_bayes_beta0,tumor_bayes_pvalue=file$tumor_bayes_pvalue)
  ftmp=file[file$unitID %in% CCid,]
  names(ftmp)=c("unitID",paste(group1[i],names(file)[2:5],sep="_"))
  rt=merge(rt,ftmp,by="unitID",all=T)
}
file=rt
j=ncol(file)-1
file$mean.normal.beta0=rowMeans(file[,seq(2,j,4)],na.rm = T)
file$mean.tumor.beta0=rowMeans(file[,seq(4,j,4)],na.rm = T)
std <- function(x) {
if(length(x[!is.na(x)])==1){return (0)
}else(return(sd(x,na.rm = T)/sqrt(length(x[!is.na(x)]))))
}###算标准误差
file$std.normal.beta0=apply(file[,seq(2,j,4)], 1, std)
file$std.tumor.beta0=apply(file[,seq(4,j,4)],1,std)
file$normal.beta0.upper=file$mean.normal.beta0+file$std.normal.beta0
file$nomal.beta0.lower=file$mean.normal.beta0-file$std.normal.beta0
file$tumor.beta0.upper=file$mean.tumor.beta0+file$std.tumor.beta0
file$tumor.beta0.lower=file$mean.tumor.beta0-file$std.tumor.beta0
#file=merge(anno,file,by="unitID")
file=merge(file,orders,by="unitID")


normaldata=select(file,unitID,mean.normal.beta0,normal.beta0.upper,nomal.beta0.lower,avsnp150,group,num)

p5=ggplot(normaldata,aes(y=mean.normal.beta0,x=num))+scale_x_continuous(breaks = normaldata$num,labels = normaldata$avsnp150,position="top")+
  geom_point(aes(color=group),size=1.5)+geom_errorbar(aes(ymin=nomal.beta0.lower,ymax=normal.beta0.upper,color=group),width=0.5)+labs(x="",y="con β0")+
  geom_hline(yintercept = 0,color="red",linetype="dashed")+scale_color_manual(values = c('#708090','red','blue'))+coord_flip()+theme_classic(base_size = 15)+
  theme(panel.grid.minor = element_blank(),panel.border = element_blank(),axis.line =element_line(size = 0.5))+theme(legend.position="none")
p5
tumordata=select(file,unitID,mean.tumor.beta0,tumor.beta0.upper,tumor.beta0.lower,avsnp150,group,num)
tumordata$na=""

p6=ggplot(tumordata,aes(y=mean.tumor.beta0,x=num))+scale_x_continuous(breaks = tumordata$num,labels = tumordata$na)+
    geom_point(aes(color=group),size=1.5)+geom_errorbar(aes(ymin=tumor.beta0.lower,ymax=tumor.beta0.upper,color=group),width=0.5)+
    labs(x="",y="case β0")+geom_hline(yintercept = 0,color="red",linetype="dashed")+scale_color_manual(values = c('#708090','red','blue'))+theme_classic(base_size = 15)+
    theme(panel.grid.minor = element_blank(),panel.border = element_blank(),axis.line = element_line(size = 0.5))+coord_flip()+theme(legend.position="none")
p6
layout <- matrix(c(1, 1, 1, 1, 2, 2, 2), nrow = 1)
multiplot(plotlist=list(p5,p6),layout=layout)





###############另开R，查看健康样本中有多少相同模式的位点
setwd("E:\\5hmc_file\\2_5hmc_yjp_bam\\ASM\\bayes_pvalue_beta0")
file=read.table("pvalue_sig_statis.csv",head=T,sep=",")
sel=which(file$num1==2&file$twins==1)
file=file[!(file$num1==2&file$twins==1),]
file=file[file$num1>1,]
id=file$unitID

group1=c("X2B_X1T","M8_M7","M6_M5","M2_M1","M48_M47","M50_M49","M28_M27","M30_M29","M26_M25","M35_M36","M18_M17","M20_M19","M22_M21","M40_M39")

i=11
fn1=paste0(group1[i],".bayes_p.txt")
file=read.table(fn1,head=T,sep = "\t")
file$unitID=paste(file$chrom,file$position,file$ref,file$var,sep=":")
file$id=paste(file$chrom,file$position,sep=":")
file=file[file$id %in% id,]
file=file[file$normal_bayes_pvalue<0.05|file$tumor_bayes_pvalue<0.05,]
file=data.frame(unitID=file$unitID,normal_bayes_beta0=file$normal_bayes_beta0,normal_bayes_pvalue=file$normal_bayes_pvalue,tumor_bayes_beta0=file$tumor_bayes_beta0,tumor_bayes_pvalue=file$tumor_bayes_pvalue)
file$normal_group="nosig"
file[file$normal_bayes_beta0>0&file$normal_bayes_pvalue<0.05,]$normal_group="up"
file[file$normal_bayes_beta0<0&file$normal_bayes_pvalue<0.05,]$normal_group="down"
file$tumor_group="nosig"
file[file$tumor_bayes_beta0>0&file$tumor_bayes_pvalue<0.05,]$tumor_group="up"
file[file$tumor_bayes_beta0<0&file$tumor_bayes_pvalue<0.05,]$tumor_group="down"
file$pattern=paste0("normal_",file$normal_group,"-tumor_",file$tumor_group)
names(file)=c("unitID",paste(group1[i],names(file)[2:5],sep="_"),names(file)[6:8])
normal_down_tumor_down=file[file$pattern=="normal_down-tumor_down",][,1:5]
one_is_down=file[file$pattern=="normal_down-tumor_nosig"|file$pattern=="normal_nosig-tumor_down",][,1:5]
onedown_oneup=file[file$pattern=="normal_down-tumor_up"|file$pattern=="normal_up-tumor_down",][,1:5]
one_is_up=file[file$pattern=="normal_nosig-tumor_up"|file$pattern=="normal_up-tumor_nosig",][,1:5]
normal_up_tumor_up=file[file$pattern=="normal_up-tumor_up",][,1:5]


for(i in 12:14){
  fn1=paste0(group1[i],".bayes_p.txt")
  file=read.table(fn1,head=T,sep = "\t")
  file$unitID=paste(file$chrom,file$position,file$ref,file$var,sep=":")
  file$id=paste(file$chrom,file$position,sep=":")
  file=file[file$id %in% id,]
  file=file[file$normal_bayes_pvalue<0.05|file$tumor_bayes_pvalue<0.05,]
  file=data.frame(unitID=file$unitID,normal_bayes_beta0=file$normal_bayes_beta0,normal_bayes_pvalue=file$normal_bayes_pvalue,tumor_bayes_beta0=file$tumor_bayes_beta0,tumor_bayes_pvalue=file$tumor_bayes_pvalue)
  file$normal_group="nosig"
  file[file$normal_bayes_beta0>0&file$normal_bayes_pvalue<0.05,]$normal_group="up"
  file[file$normal_bayes_beta0<0&file$normal_bayes_pvalue<0.05,]$normal_group="down"
  file$tumor_group="nosig"
  file[file$tumor_bayes_beta0>0&file$tumor_bayes_pvalue<0.05,]$tumor_group="up"
  file[file$tumor_bayes_beta0<0&file$tumor_bayes_pvalue<0.05,]$tumor_group="down"
  file$pattern=paste0("normal_",file$normal_group,"-tumor_",file$tumor_group)
  names(file)=c("unitID",paste(group1[i],names(file)[2:5],sep="_"),names(file)[6:8])
  normal_down_tumor_down=merge(normal_down_tumor_down,file[file$pattern=="normal_down-tumor_down",][,1:5],by="unitID")
  one_is_down=merge(one_is_down,file[file$pattern=="normal_down-tumor_nosig"|file$pattern=="normal_nosig-tumor_down",][,1:5],by="unitID")
  onedown_oneup=merge(onedown_oneup,file[file$pattern=="normal_down-tumor_up"|file$pattern=="normal_up-tumor_down",][,1:5],by="unitID")
  one_is_up=merge(one_is_up,file[file$pattern=="normal_nosig-tumor_up"|file$pattern=="normal_up-tumor_nosig",][,1:5],by="unitID")
  normal_up_tumor_up=merge(normal_up_tumor_up,file[file$pattern=="normal_up-tumor_up",][,1:5],by="unitID")
}

library(dplyr)
library(ggplot2)
normal_up_tumor_up$group="normal_up_tumor_up"
normal_down_tumor_down$group="normal_down_tumor_down"
one_is_down$group="one_is_down"
file=rbind(normal_up_tumor_up,normal_down_tumor_down,one_is_down)
j=ncol(file)-1
file$mean.normal.beta0=rowMeans(file[,seq(2,j,4)])
file$mean.tumor.beta0=rowMeans(file[,seq(4,j,4)])
std <- function(x) sd(x)/sqrt(length(x))###算标准误差
file$std.normal.beta0=apply(file[,seq(2,j,4)], 1, std)
file$std.tumor.beta0=apply(file[,seq(4,j,4)],1,std)
file$normal.beta0.upper=file$mean.normal.beta0+file$std.normal.beta0
file$nomal.beta0.lower=file$mean.normal.beta0-file$std.normal.beta0
file$tumor.beta0.upper=file$mean.tumor.beta0+file$std.tumor.beta0
file$tumor.beta0.lower=file$mean.tumor.beta0-file$std.tumor.beta0
anno=read.csv("../20201107/HC.intersect.anno.hg19_multianno.csv",head=T)
anno$unitID=paste(anno$Chr,anno$Start,anno$Ref,anno$Alt,sep=":")
file=merge(anno,file,by="unitID")

normaldata=select(file,unitID,mean.normal.beta0,normal.beta0.upper,nomal.beta0.lower,avsnp150,group)
normaldata=arrange(normaldata,group,mean.normal.beta0)
normaldata$num=1:dim(normaldata)[1]
orders=select(normaldata,unitID,num)

p1=ggplot(normaldata,aes(y=mean.normal.beta0,x=num))+scale_x_continuous(breaks = normaldata$num,labels = normaldata$avsnp150,position="top")+
  geom_point(aes(color=group),size=1.5)+geom_errorbar(aes(ymin=nomal.beta0.lower,ymax=normal.beta0.upper,color=group),width=0.5)+labs(x="",y="con β0")+
  geom_hline(yintercept = 0,color="red",linetype="dashed")+scale_color_manual(values = c('#708090','red','blue'))+coord_flip()+theme_classic(base_size = 15)+
  theme(panel.grid.minor = element_blank(),panel.border = element_blank(),axis.line =element_line(size = 0.5))+theme(legend.position="none")
p1
tumordata=select(file,unitID,mean.tumor.beta0,tumor.beta0.upper,tumor.beta0.lower,avsnp150,group)
tumordata$na=""
tumordata=merge(tumordata,orders,by="unitID")

p2=ggplot(tumordata,aes(y=mean.tumor.beta0,x=num))+scale_x_continuous(breaks = tumordata$num,labels = tumordata$na)+
  geom_point(aes(color=group),size=1.5)+geom_errorbar(aes(ymin=tumor.beta0.lower,ymax=tumor.beta0.upper,color=group),width=0.5)+
  labs(x="",y="case β0")+geom_hline(yintercept = 0,color="red",linetype="dashed")+scale_color_manual(values = c('#708090','red','blue'))+theme_classic(base_size = 15)+
  theme(panel.grid.minor = element_blank(),panel.border = element_blank(),axis.line = element_line(size = 0.5))+coord_flip()+theme(legend.position="none")
p2
layout <- matrix(c(1, 1, 1, 1, 2, 2, 2), nrow = 1)
multiplot(plotlist=list(p1,p2),layout=layout)

HCid=file$unitID
original.normaldata=normaldata
orders=data.frame(unitID=original.normaldata$unitID,original.normaldata[,5:7])

###看单发样本中的
i=1
fn1=paste0(group1[i],".bayes_p.txt")
file=read.table(fn1,head=T,sep = "\t")
file$unitID=paste(file$chrom,file$position,file$ref,file$var,sep=":")
file=data.frame(unitID=file$unitID,normal_bayes_beta0=file$normal_bayes_beta0,normal_bayes_pvalue=file$normal_bayes_pvalue,tumor_bayes_beta0=file$tumor_bayes_beta0,tumor_bayes_pvalue=file$tumor_bayes_pvalue)
rt=file[file$unitID %in% HCid,]
names(rt)=c("unitID",paste(group1[i],names(file)[2:5],sep="_"))

for (i in 2:6) {
  fn1=paste0(group1[i],".bayes_p.txt")
  file=read.table(fn1,head=T,sep = "\t")
  file$unitID=paste(file$chrom,file$position,file$ref,file$var,sep=":")
  file=data.frame(unitID=file$unitID,normal_bayes_beta0=file$normal_bayes_beta0,normal_bayes_pvalue=file$normal_bayes_pvalue,tumor_bayes_beta0=file$tumor_bayes_beta0,tumor_bayes_pvalue=file$tumor_bayes_pvalue)
  ftmp=file[file$unitID %in% HCid,]
  names(ftmp)=c("unitID",paste(group1[i],names(file)[2:5],sep="_"))
  rt=merge(rt,ftmp,by="unitID",all=T)
}
file=rt
j=ncol(file)-1
file$mean.normal.beta0=rowMeans(file[,seq(2,j,4)],na.rm = T)
file$mean.tumor.beta0=rowMeans(file[,seq(4,j,4)],na.rm = T)
std <- function(x) {
if(length(x[!is.na(x)])==1){return (0)
}else(return(sd(x,na.rm = T)/sqrt(length(x[!is.na(x)]))))
}###算标准误差
file$std.normal.beta0=apply(file[,seq(2,j,4)], 1, std)
file$std.tumor.beta0=apply(file[,seq(4,j,4)],1,std)
file$normal.beta0.upper=file$mean.normal.beta0+file$std.normal.beta0
file$nomal.beta0.lower=file$mean.normal.beta0-file$std.normal.beta0
file$tumor.beta0.upper=file$mean.tumor.beta0+file$std.tumor.beta0
file$tumor.beta0.lower=file$mean.tumor.beta0-file$std.tumor.beta0
#file=merge(anno,file,by="unitID")
file=merge(file,orders,by="unitID")


normaldata=select(file,unitID,mean.normal.beta0,normal.beta0.upper,nomal.beta0.lower,avsnp150,group,num)

p3=ggplot(normaldata,aes(y=mean.normal.beta0,x=num))+scale_x_continuous(breaks = normaldata$num,labels = normaldata$avsnp150,position="top")+
  geom_point(aes(color=group),size=1.5)+geom_errorbar(aes(ymin=nomal.beta0.lower,ymax=normal.beta0.upper,color=group),width=0.5)+labs(x="",y="con β0")+
  geom_hline(yintercept = 0,color="red",linetype="dashed")+scale_color_manual(values = c('#708090','red','blue'))+coord_flip()+theme_classic(base_size = 15)+
  theme(panel.grid.minor = element_blank(),panel.border = element_blank(),axis.line =element_line(size = 0.5))+theme(legend.position="none")
p3
tumordata=select(file,unitID,mean.tumor.beta0,tumor.beta0.upper,tumor.beta0.lower,avsnp150,group,num)
tumordata$na=""

p4=ggplot(tumordata,aes(y=mean.tumor.beta0,x=num))+scale_x_continuous(breaks = tumordata$num,labels = tumordata$na)+
    geom_point(aes(color=group),size=1.5)+geom_errorbar(aes(ymin=tumor.beta0.lower,ymax=tumor.beta0.upper,color=group),width=0.5)+
    labs(x="",y="case β0")+geom_hline(yintercept = 0,color="red",linetype="dashed")+scale_color_manual(values = c('#708090','red','blue'))+theme_classic(base_size = 15)+
    theme(panel.grid.minor = element_blank(),panel.border = element_blank(),axis.line = element_line(size = 0.5))+coord_flip()+theme(legend.position="none")
p4
layout <- matrix(c(1, 1, 1, 1, 2, 2, 2), nrow = 1)
multiplot(plotlist=list(p3,p4),layout=layout)

###看双发样本中的
i=7
fn1=paste0(group1[i],".bayes_p.txt")
file=read.table(fn1,head=T,sep = "\t")
file$unitID=paste(file$chrom,file$position,file$ref,file$var,sep=":")
file=data.frame(unitID=file$unitID,normal_bayes_beta0=file$normal_bayes_beta0,normal_bayes_pvalue=file$normal_bayes_pvalue,tumor_bayes_beta0=file$tumor_bayes_beta0,tumor_bayes_pvalue=file$tumor_bayes_pvalue)
rt=file[file$unitID %in% HCid,]
names(rt)=c("unitID",paste(group1[i],names(file)[2:5],sep="_"))

for (i in 8:10) {
  fn1=paste0(group1[i],".bayes_p.txt")
  file=read.table(fn1,head=T,sep = "\t")
  file$unitID=paste(file$chrom,file$position,file$ref,file$var,sep=":")
  file=data.frame(unitID=file$unitID,normal_bayes_beta0=file$normal_bayes_beta0,normal_bayes_pvalue=file$normal_bayes_pvalue,tumor_bayes_beta0=file$tumor_bayes_beta0,tumor_bayes_pvalue=file$tumor_bayes_pvalue)
  ftmp=file[file$unitID %in% HCid,]
  names(ftmp)=c("unitID",paste(group1[i],names(file)[2:5],sep="_"))
  rt=merge(rt,ftmp,by="unitID",all=T)
}
file=rt
j=ncol(file)-1
file$mean.normal.beta0=rowMeans(file[,seq(2,j,4)],na.rm = T)
file$mean.tumor.beta0=rowMeans(file[,seq(4,j,4)],na.rm = T)
std <- function(x) {
if(length(x[!is.na(x)])==1){return (0)
}else(return(sd(x,na.rm = T)/sqrt(length(x[!is.na(x)]))))
}###算标准误差
file$std.normal.beta0=apply(file[,seq(2,j,4)], 1, std)
file$std.tumor.beta0=apply(file[,seq(4,j,4)],1,std)
file$normal.beta0.upper=file$mean.normal.beta0+file$std.normal.beta0
file$nomal.beta0.lower=file$mean.normal.beta0-file$std.normal.beta0
file$tumor.beta0.upper=file$mean.tumor.beta0+file$std.tumor.beta0
file$tumor.beta0.lower=file$mean.tumor.beta0-file$std.tumor.beta0
#file=merge(anno,file,by="unitID")
file=merge(file,orders,by="unitID")


normaldata=select(file,unitID,mean.normal.beta0,normal.beta0.upper,nomal.beta0.lower,avsnp150,group,num)

p5=ggplot(normaldata,aes(y=mean.normal.beta0,x=num))+scale_x_continuous(breaks = normaldata$num,labels = normaldata$avsnp150,position="top")+
  geom_point(aes(color=group),size=1.5)+geom_errorbar(aes(ymin=nomal.beta0.lower,ymax=normal.beta0.upper,color=group),width=0.5)+labs(x="",y="con β0")+
  geom_hline(yintercept = 0,color="red",linetype="dashed")+scale_color_manual(values = c('#708090','red','blue'))+coord_flip()+theme_classic(base_size = 15)+
  theme(panel.grid.minor = element_blank(),panel.border = element_blank(),axis.line =element_line(size = 0.5))+theme(legend.position="none")
p5
tumordata=select(file,unitID,mean.tumor.beta0,tumor.beta0.upper,tumor.beta0.lower,avsnp150,group,num)
tumordata$na=""

p6=ggplot(tumordata,aes(y=mean.tumor.beta0,x=num))+scale_x_continuous(breaks = tumordata$num,labels = tumordata$na)+
    geom_point(aes(color=group),size=1.5)+geom_errorbar(aes(ymin=tumor.beta0.lower,ymax=tumor.beta0.upper,color=group),width=0.5)+
    labs(x="",y="case β0")+geom_hline(yintercept = 0,color="red",linetype="dashed")+scale_color_manual(values = c('#708090','red','blue'))+theme_classic(base_size = 15)+
    theme(panel.grid.minor = element_blank(),panel.border = element_blank(),axis.line = element_line(size = 0.5))+coord_flip()+theme(legend.position="none")
p6
layout <- matrix(c(1, 1, 1, 1, 2, 2, 2), nrow = 1)
multiplot(plotlist=list(p5,p6),layout=layout)

